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1A»  ABSTRACT 

Level  Set  Systems  (LSS)  has  developed  a  3D  Euteflan,  time  domain  method  forcompidlng  the  propagation  and  scattering  of  short 
electromagnetic  pulses  over  arbitrarily  tong  distances.  Efforts  accommodated  Include  variation  of  Indexof  refraction  and  scattering 
from  complex  featires  including  tanks  and  aircraft.  LSS  developed  a  level  set  method  for  the  construction  of  wavefronts  that  handles 

resolution,  multivalued  solutions,  reflection  and  refraction.  The  key  idea  involves 
an  Bulerian,  fixed  grid  solver  involving  a  system  of  linear  Liouville  equations,  three 
equations  in  five  space  dimensions .  The  complexity  of  each  update  is  no  worse  than  that 
of  ray  tracing  because  of  the  use  of  a  new  local  level  set  method  for  high  dimension 
and  codimension.  This  allows  this  level  set  method  to  be  a  viable  alternative  to  ray 
tracing  with  many  advantages  for  realistic  DOD  applications  involving  geometric  optics. 
Flow  Analysis  Inc  (FAI)  has  developed  a  lattice  confinement  method  that  has  been  shown 
To  be  effective  for  propagating  short  pulses  as  nonlinear  solitary  waves  over  long 
distances.  These  two  methods  complement  each  other  and  lay  the  foundation  for  greatly 
improved  Radar  Cross  Section  computations. _ 
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1.  INTRODUCTION 

Geometric  optics  makes  its  impact,  both  in  mathematics  and  in  real  world 
applications  related  to  ray  tracing,  migration  and  tomography.  The  primary  objective  of 
this  phase  of  the  research  has  been  to  develop  a  3-D  Eulerian  time  domain  method  for 
accurately  computing  the  propagation  and  scattering  of  short  electromagnetic  and 
acoustic  waves  over  very  long  distances  ( » 104  A )  through  the  atmosphere,  where  A,  is 
the  width  of  the  pulse.  Effects  accommodated  include  variation  of  index  of  refraction  and 
scattering  from  complex  features,  such  as  tanks  and  aircraft.  Of  special  interest  in  these 
problems  are  the  wavefronts,  or  points  of  constant  travel  time  away  from  sources,  in  the 
medium.  Previously  in  [OCKST]  under  the  support  of  this  contract,  we  initiated  a  level 
set  approach  for  the  construction  of  wavefronts  in  isotropic  media  that  handled  the  two 
major  algorithmic  issues  of  resolution  and  multivalued  solutions.  This  approach  was  quite 
general  and  we  were  able  to  construct  wavefronts  in  the  presence  of  refraction,  reflection, 
higher  dimensions  and,  in  other  new  work  supported  by  this  contract  [QCO,  COQ],  we 
included  anisotropy. 

However,  the  technique  proposed  for  handling  reflections  of  waves  off  objects,  the 
basic  phenomenon  involved  in  all  applications  of  geometric  optics,  was  inefficient  and 
unwieldy  to  the  point  of  being  unusable,  especially  in  the  presence  of  multiple 
reflections.  In  [CKOST]  (enclosed)  supported  by  this  contract  we  introduced  an  approach 
based  on  the  foundation  presented  in  [OCKST]  that  fixes  this  issue.  This  reworking  fully 
allows  this  level  set  method  to  be  a  viable  alternative  to  ray  tracing,  with  many 
advantages,  for  realistic  applications  involving  geometric  optics,  especially  applications 
of  interest  to  DOD. 

The  main  difficulty  encountered  by  all  numerical  approaches  in  the  construction  of 
wavefronts  in  the  geometric  optics  or  ray  tracing  setting  lies  in  a  choice  of  either 
resolution  or  generation  of  multivalued  solutions.  Multivaluedness  in  wavefronts  occur 
when  wavefronts  cross  and  more  than  one  ray  occupies  a  region  in  space.  The  formation 
of  the  well  known  phenomenon  of  caustics  originates  from  this.  Solutions  obtained 
following  the  Lagrangian  representation  of  ray  tracing,  which  involves  using  the  method 
of  characteristics  to  track  the  position  and  ray  direction  of  points  on  the  wavefronts,  can 
automatically  produce  multivalued  wavefronts  but  have  difficulties  in  resolving 
wavefronts  in  general,  especially  when  they  are  diverging.  On  the  other  hand,  Eulerian 
approaches  applied  to  the  eikonal  equation  automatically  resolve  all  wavefronts  over  an 
underlying  grid  in  space  but  encounter  difficulties  in  generating  multivalued  wavefronts. 
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We  overcame  these  difficulties  in  our  first  paper  [OCKST],  and  now  have  extended  this 
to  anisotropy  and  reflection.  In  [JO],  we  extended  the  idea  to  the  computation  of 
multivalued  solutions  to  quasi-linear  hyperbolic  partial  equations  and  Hamilton- Jacobi 
equations  in  arbitrary  space  dimensions.  We  use  the  classic  idea  of  Courant  and  Hilbert  to 
define  the  solution  of  the  quasi-linear  hyperbolic  PDEs  or  the  gradient  of  the  solution  to 
the  Hamilton- Jacobi  equations  as  zero  level  sets  of  level  set  functions.  Then  the  evolution 
equations  for  the  level  set  functions  satisfy  linear  Liouville  equations  defined  in  the 
“phase”  space,  unfolding  the  singularities  and  preventing  the  numerical  capturing  of 
viscosity  solution.  This  provides  a  computational  framework  for  the  computations  of 
multivalued  geometric  solutions  to  general  quasilinear  PDEs.  In  [CLO],  we  apply  our 
idea  to  compute  the  linear  Schrodinger  equations  with  efficiently  highly  oscillating  initial 
data.  The  high-frequency  asymptotics  of  these  equations  lead  to  the  well  known  WKB 
system,  where  the  phase  evolves  according  to  a  nonlinear  Hamilton- Jacobi  equation  and 
the  intensity  of  the  plane  wave  is  governed  by  a  linear  conservation  law.  This  becomes 
equivalent  to  solving  a  linear  Liouville  equation  with  multi-valued  phase  as  we  already 
developed  in  [OCKST]. 

In  higher  dimensions,  computational  speed,  visualization  and  efficiency  start  to 
become  major  issues.  We  have  begun  successfully  attacking  these  issues  and  have  a  local 
method,  involving  local  memory  and  an  innovative  new  visualization  technique  for 
submanifolds  of  high  dimensional  space  developed  by  C.  Min  and  S.  Osher[M,  MO]. 
Also,  we  are  investigating  a  suggestion  put  forward  by  F.  Reitich  and  J.  Qian  that  spectral 
methods  be  used  to  update  our  evolution  in  the  angle  variables  and  discontinuous 
Galerkin  methods  be  used  to  update  the  space  variables.  Qian  and  Reitich  are  using  our 
basic  technique  and  we  expect  to  work  with  them. 

Another  advantage  of  our  Eulerian  grid  based  techniques  is  that  it  can  be  updated 
simultaneously  with  full  Maxwell’s  equation  solvers.  We  intend  to  work  with  our 
colleagues  at  Hypercomp,  Inc.  and  Brown  University  on  this.  Finally,  we  are  exploring  a 
promising  new  confinement  technique,  put  forward  by  our  subcontractor,  Flow  Analysis 
Inc  to  help  enhance  both  of  these  methods. 

The  Flow  Analysis,  Inc.’s  pulse  propagation  algorithm  involves  a  new  -  “Lattice 
Confinement”  Method  (LCM)  that  has  been  shown  to  be  effective  for  propagating  short 
pulses  ( ~  2  cells  wide)  as  non  linear  solitary  waves  over  arbitrary  distances  with  no 
numerical  spreading.  This  pulse  algorithm,  by  itself,  will  be  an  important  contribution  to 
time  domain  radar  scattering  techniques. 

The  short  pulses  resulting  from  the  LCM  have  been  previously  shown  to  also 
allow  the  computation  of  high  frequency  time  domain  phases  and  amplitudes,  including 
multiple  arrival  times.  This  aspect  represents  an  important  capability.  In  this  mode, 
effectively,  as  the  short  pulse  surface  sweeps  through  the  computational  grid,  new 
amplitude/phase  values  will  be  generated  within  the  band  extending  across  the  width  of 
the  pulse.  After  passage  of  the  pulse,  these  values  will  be  stored  in  an  array  so  that  new 
pulses,  and  hence  new  arriving  waves  can  be  treated.  This  capability  has  previously  been 
presented  in  [SWUP]. 
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1.1  Relation  to  Other  Methods 


a)  Direct  Wave  Equation  Solvers 

The  LCM  will  be  useful  by  itself,  but  will  also  be  able  to  be  matched  to  other, 
existing  methods/codes,  such  as  direct  Maxwell  equation  solvers  [S],  which  are  in  use  for 
computing  the  fields  near  an  object.  In  this  mode,  the  LCM  will  serve  to  extend  the 
solvers  far  beyond  their  limit.  Since  direct  solvers  resolve  the  individual  waves  (even 
optimized  methods  [CSE]  still  must  satisfy  the  Nyquist  criterion),  the  computational 
requirements  limit  them  to  regions  of  0(1 02 A).  Thus  these  methods,  alone,  are  not 
sufficient  to  treat  long  distance  propagation,  either  in  the  frequency  or  time  domain.  This 
matching  feature  should  be  important  for  low  observable  objects  where  complex,  well- 
developed  computational  methods  are  already  in  use  near  the  objects. 

(b)  Eikonal  Schemes 

As  described  above,  the  LCM  will  be  able  to  automatically  transition  to 
conventional  Eulerian  solves  near  a  scatterer,  if  necessary.  This  appears  to  be  difficult 
with  conventional  high  frequency  Eikonal  methods.  Also,  as  described  in  [SWUP], 
multiple-arrival  times  will  automatically  be  accommodated,  something  difficult  for 
Eikonal  schemes  [FEO],  but  which  is  easy  for  the  vector  level  set  method  of  LSS. 

In  this  way,  the  LCM  complements  the  other  research  discussed  in  this  report. 

(c)  Kirchhoff  Methods 

Variable  index  of  refraction  effects  will  be  treated,  which  are  not  included  in 
conventional  Kirchhoff  integral  methods. 

(d)  Lippmann  Schwinger  Integral  Schemes 

Almost  all  of  the  problems  to  be  treated  involve  waves  that  pass  through  each 
point  either  once,  or  only  a  few  times.  This  includes  multiple  scattering  from  surfaces  and 
refraction/ducting  effects.  In  many  of  these,  the  physical  space  can  be  divided  into 
separate  regions  that  can  be  computed  sequentially  -  such  as  the  region  near  a  scatter  and 
the  far  field.  This  feature  can  be  used  by  the  LCM  to  easily  decompose  the  computational 
domain  into  separate  domains  to  be  treated  sequentially,  thereby  greatly  reducing 
memory  requirements  and  allowing  parallel  computation.  This  does  not  appear  to  be 
feasible  with  straightforward  implementations  of  the  Lippmann  Schwinger  integral 
equation  approaches  [BK],  where  each  region  of  space  is  connected  to  other  regions  by  a 
Green’s  function  (assuming  varying  index  of  refraction).  Further,  with  the  LCM, 
computations  are  done  only  in  a  narrow  band  of  the  grid  ( n  ~  0(1)  cells)  at  each  of  a 
number  of  time  steps  (NT).  This  number  is  of  the  order  of  the  grid  size  in  each  direction 
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(say  Nx ).  Thus,  for  3-D  computations,  the  total  number  of  computations  is  only  0{N\ ) . 
Also,  each  of  these  computations  is  a  simple,  local  finite  difference  operation. 

(e)  General  Frequency  Domain  Methods 

The  LCM  can  be  used  to  compute  time  domain  solutions  for  long-range 
propagation/scattering  of  short  pulses,  which  can  be  important  for  modem  radar.  This 
appears  to  be  outside  the  scope  of  narrow  band  frequency  domain  schemes,  since  short 
pulses  consist  of  a  large  frequency  band  which  would  require  many  frequency  domain 
solutions  using  conventional  methods. 

(f)  Ray  Tracing  Methods 

Ray  tracing  schemes  have  well-known  difficulties  in  reconstructing  continuous 
pulse  or  phase  surfaces  in  general.  Also,  information  must  be  interpolated  onto  a  uniform 
grid  to  be  useful  for  many  purposes. 

The  LCM,  as  are  methods  (a)  -  (e)  above,  is  completely  Eulerian.  A  convenient 
structured  computational  grid  is  used  which  does  not  have  any  special  data  handling 
requirements.  This  grid  is  suitable  for  representing  variations  of  the  index  of  refraction, 
topographic  boundary  conditions  and  the  geometries  of  scattering  objects.  Also, 
transmission  and  scattering  amplitudes  and  phases,  if  they  are  to  be  useful,  should  be 
represented  on  a  convenient,  structured  grid. 

The  following  Level  Set  Systems,  Inc.  personnel  were  supported  under  this  grant: 
Dr.  Stanley  Osher,  Dr.  Susan  Chen,  Dr.  Hyeseon  Shim  and  Dr.  Myungjoo  Kang.  The 
following  Flow  Analysis,  Inc.  personnel  were  supported  under  this  grant:  Dr.  John 
Steinhoff  and  Dr.  William  Dietz. 


2.  LEVEL  SET  FORMULATION  AND  VISUALIZATION 


For  purpose  of  simplicity  and  exposition,  we  restrict  to  our  attention  to  n  =  2  for 
formulation,  but  we  will  give  the  result  for  higher  dimension  too.  Thus  reduced  phase 
space,  for  a  fixed  time,  can  be  written  as  the  set  of  (x,  9) ,  where  x  e  R2  and  6  e  \-n,  n] , 


9  being  the  angle  of  p  in  polar  coordinates.  So  for  a  fixed  time,  the  representation  of  the 
wavefront,  called  bicharacteristic  strip,  is  a  smooth  curve  in  reduced  phase  space.  In  our 
formulation,  the  curve  is  represented  by  the  intersection  of  two  surfaces,  which  in  turn  are 
represented  by  the  zero  level  sets  of  two  level  set  functions.  Denoting  the  level  set 
functions  by  </>  and  iff ,  the  curve  is  thus  the  set  of  points  where  <f>  =  if/  =  0 .  The  Liouville 
equation  on  each  level  set  function,  written  as 

^+V*v*^  =  ° 

V'/+V*v,^  =  ° 


(2.1) 


where  9  denotes  the  phase  plane  angle  and  v  is  given  by 
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f 


v(x,  9)  = 


c(x)cos# 

c(x)sin# 


(2.2) 


^  (x)  sin  9  -  cXi  (x)  cos  9  J 

where  c  is  the  given  local  wave  velocity  permitted  in  the  medium.  Also  note  these  two 
transport  equations  in  <j)  and  if/  can  be  solved  separately.  For  the  case  of  refraction,  we 
use  the  smooth  delta  function  between  two  medium  to  avoid  the  sudden  jump  which  may 
create  the  instability  in  numerical  calculation.  For  the  case  of  reflection,  we  consider  it  as 
an  initial  boundary  value  problem,  i.e.,  we  give  the  boundary  condition  on  the  reflected 
wall  whenever  it  needed  and  solve  the  Liouville  equation  with  those  boundary  conditions. 
In  3  dimensional  cases,  we  can  still  use  a  uniform  grid  without  any  problem.  But  for 
higher  dimensional  cases  (n  >  2),  it  is  almost  impossible  to  use  a  uniform  mesh  to 
simulate  our  optics  problems.  The  reason  is  that  we  need  humongous  memory  and  time  to 
get  a  piece  of  result  in  higher  dimension.  One  more  issue  in  higher  dimension  cases  is 
how  to  visualize  our  results.  For  n  =  3  (i.e.,  we  need  5  dimensional  calculation),  our  result 
surface  is  the  intersection  of  three  4  dimensional  zero  level  isosurfaces. 

In  order  to  overcome  those  issues,  we  used  a  new  local  level  set  method  and 
visualization  technique  developed  by  C.  Min  and  S.  Osher  [M,  MO].  We  can  reduce  the 


numerical  cost  for  both  memory  and  time  as  much  as  0(n2)  ratio  for  n  =  3. 


3.  The  Lattice  Confinement  Method 


3.1  Basic  Features 

The  main  requirement  of  the  method  is  to  compute  thin  propagating  wave 
equation  pulses  that  do  not  spread  due  to  numerical  effects. 

The  LCM  efficiently  simulates  these  wave  equation  pulses,  in  Eulerian 
computations  on  fixed,  coarse  grids.  The  method  involves  treating  the  features  as  solitary 
waves  that  obey  nonlinear,  difference  equations,  which  are  different  from  Taylor 
expansion-  based  discrete  approximations  to  the  governing,  partial  differential  equations 
(pde’s).  These  equations  are  rotationally  invariant  generalizations  to  multiple  dimensions 
of  1-D  discontinuity  confinement  schemes.  The  method  is  a  generalization  of  an  earlier 
method,  “Vorticity  Confinement”  [SU],  which  was  successful  in  efficiently  treating  thin, 
vortical  regions. 

For  long  distance  propagation  of  pulses,  direct  discretization  and  solution  of  the 
governing  partial  differential  equations  (pde’s)  using  conventional  Eulerian  Taylor 
expansion-based  numerical  methods  to  resolve  the  thin  features  can  be  prohibitively 
expensive.  Even  adaptive  unstructured  grid  methods  are  very  expensive  and  complex  for 
general  problems  with  many  small-scale,  time-dependent  features.  Fortunately,  the  details 
of  the  internal  structure  of  these  small  features  are  often  not  as  important  as  integral 
quantities.  The  quantities  of  importance  for  our  purpose  are  the  centroid  motion  and  total 
integrated  amplitude  at  each  point  along  the  pulse  surface.  The  main  issue  in  computing 
these  cases  is  that  conventional  pde-based  methods  require  a  relatively  large  number  of 
grid  cells  (4-8)  across  each  small  dimension  to  treat  a  feature,  such  as  a  wave  equation 
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pulse.  Even  then,  the  details  of  the  internal  structure  will  be  mainly  determined  by  the 
discrete  numerics,  and  not  the  physics  of  the  pde.  Also,  numerical  discretization  errors 
will  still  build  up  over  long  distances,  causing,  for  example,  large,  unphysical  spreading. 

This  leads  us  to  the  idea  of  simulating  or  “modeling  “  the  thin  features  directly  on 
the  grid  with  difference  equations,  rather  than  attempting  to  accurately  discretize  the 
pde’s  for  them  using  finite  difference  approximations.  The  idea  of  modeling,  or  solving 
for  small-scale  features  directly  on  the  grid  without  using  smoothness  assumptions  or 
Taylor  expansions,  i.e.,  as  “weak  solutions”,  goes  back  to  work  of  Lax  and  others  [L],  but 
was  applied  mostly  to  shocks.  Shocks,  however,  effectively,  “capture  themselves” 
because  they  have  converging  characteristics.  Harten  [H]  did  treat  contact  discontinuities 
in  this  way,  which  do  not  have  converging  characteristics,  but  for  1-D  compressible  flow. 

A  major  difference  between  the  LCM  and  conventional  discontinuity  steepening 
schemes  is  that  LCM  is  not  a  concatenation  of  1-D  operators  along  each  axis,  but  is 
intrinsically  rotationally  invariant. 

First,  the  Lattice  Confinement  Method  will  be  reviewed  for  wave  equation  pulses. 
Some  results  will  then  be  presented  for  short  scalar  pulses.  Previous  results  [SDHXLF] 
for  a  2-D  pulse  reflecting  from  planar  surfaces  will  be  shown,  then  3-D  results  for  a  pulse 
reflecting  from  complex  objects  (a  missile  and  an  aircraft).  The  missile  case  will  involve 
multiple,  curvilinear  grids. 

3.2  Basic  LCM  Algorithm 

(Parts  of  this  section  are  related  to  Ref.  [SDHXLF]) 

As  mentioned,  the  main  idea  is  to  treat  thin  features  as  nonlinear  solitary  waves 
that  “live”  directly  on  the  grid  lattice,  spread  over  only  a  few  cells.  The  internal  structure 
is  determined  by  the  discrete  lattice  equations.  The  total  amplitude,  centroid  and  (in  a 
future  extension),  a  few  moments,  however,  are  determined  by  the  physics.  These 
quantities  are  transported  across  the  grid  with  essentially  no  numerical  errors. 

Essentially,  these  discrete  equations  define  a  simple,  implicit  model  which  obeys 
a  “fast”  dynamics,  relaxing  to  an  asymptotic,  propagating  state  in  a  few  time  steps.  In  this 
way  we  can  simulate  the  most  important  physical  effects  of  the  small  scales,  which 
cannot  be  accurately  computed  by  just  discretizing  the  governing  pde’s  on  the  given  grid. 
These  effects  can  be,  for  example,  that  thin  wave  equation  pulses  propagate  over  long 
distances  without  spreading  in  a  smooth,  slowly  varying  external  field,  and  that  they  can 
merge  or  reflect,  respectively,  and  thus  change  topology. 

3.2.1  Stationary  Case 

The  formulation  presented  here  is  related  to  that  presented  in  [FWS]  in  1-D.  First 
a  stationary  pulse  is  discussed,  requiring  only  an  iteration  of  the  Lattice  Confinement 
terms,  so  that  the  simple  asymptotic  form  can  be  seen.  Advection,  in  general,  will  change 
this  form  somewhat.  However,  in  the  limit  of  small  advection  time  step,  or  if  a  number  of 
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these  “Confinement”  steps  are  taken  for  each  convection  step,  then  the  following  form 
should  result.  The  same  is  true  for  the  wave  equations.  Results  very  close  to  these  are  also 
found  with  advection  steps  that  are  not  small;  these  are  shown  in  Ref.  [FWS]. 

We  start  with  an  iteration  for  a  single-signed  scalar,  </> : 

=  “V2  -  £d)] 

At 

or 

f+'  =f  +h2V2(^n -e <D")  (3.1) 

where  O  is  a  nonlinear  function  of  <f>  (given  below)  and  fj.  is  a  diffusion  coefficient  that 
can  include  numerical  effects  in  a  convection  or  wave  equation  solution  (we  assume 
physical  diffusion  is  much  smaller).  The  discretized  grid  cell  size  is  h  and  time  step,  At. 
For  the  last  term,  s  is  a  numerical  coefficient  that,  together  with  fi ,  controls  the  size  and 
time  scales  of  the  confined  features.  For  this  reason,  we  refer  to  the  two  terms  in  the 
brackets  as  “Confinement  terms”. 

The  two  (positive)  parameters,  e  and  /u ,  are  determined  by  the  two  small  scales 

of  the  computation,  h  and  At,  since  we  want  the  small  features  to  relax  to  their  solitary 
wave  shape  in  a  small  number  of  time  steps  and  to  have  a  support  of  a  small  number  of 
grid  cells.  Thus,  even  though  h  may  be  small,  the  Laplacian  will  be  large  and  the  total 
effect  also  large. 

For  the  propagating  pulse  problem,  it  is  assumed  that  the  slowness  field  is  slowly 
varying  compared  to  these  scales.  Relaxation  of  this  assumption  to  accommodate 
discontinuity  is  straightforward.  We  then  have  a  two-scale  problem  with  the  thin  structure 
obeying  a  “fast”  dynamics: 

fO  «  0 . 

With  propagation  in  a  “slow”,  smooth  external  field,  this  relation  is  still  approximately 
satisfied,  as  verified  by  computations  and  heuristic  arguments  [SDHXLF]. 

There  are  many  possibilities  for  O .  A  simple  class  is 

£c,  (#*)-' 

'  sc, 

/ 

$n  =1  <t>  I”  +s . 

The  above  sum  is  over  a  set  of  grid  nodes  near  and  including  the  node  where  O  is 
computed.  The  absolute  value  is  taken  and  S ,  a  small  positive  constant  (~  10-8 ),  is  added 
to  prevent  problems  due  to  finite  precision.  The  coefficients,  C, ,  can  depend  on  /,  but 
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good  results  for  many  cases  are  obtained  by  simply  setting  them  as  well  as  p  to  1 .  Then, 
<t>  is  the  harmonic  mean  of  <j)  on  the  local  stencil.  Other  forms  could  also  be  used,  with 
p>  1.  p  =  oo  corresponds  to  the  minimum  of  the  absolute  value:  for  2-D  and  3-D 
applications  discontinuous  operators  such  as  “min”  will  not  result  in  as  smooth 
distributions  as  continuous  ones,  and  we  use  only  p  =  1  or  p  =  2 . 

An  important  feature  of  Lattice  Confinement  is  that  all  terms  are  homogeneous  of 
degree  1  in  Eqn.  3.1  (as  they  are  in  the  original  wave  equation).  This  is  necessary  because 
the  Lattice  Confinement  terms  should  not  depend  on  the  scale  of  the  quantity  being 
confined.  Another  important  point  is  that  wavelengths  longer  than  the  thin  features  that 
are  to  be  confined  must  have  a  negative  diffusive  behavior,  so  that  the  features  remain 
confined,  that  is  stable  to  perturbations  against  spreading.  This  means  that  <D  must  be 
nonlinear:  It  is  easy  to  show  by  Von  Neumann  analysis  that  a  linear  combination  of 
terms,  for  example  of  second  and  fourth  order,  cannot  lead  to  a  stable  Confinement  for 
any  finite  range  of  coefficients:  any  wavelength  that  exhibits  negative  diffusion  would 
then  eventually  diverge. 

3.3  Wave  Equation  Formulation 

We  start  with  the  1-D  scalar  wave  equation  with  constant  wave  speed,  c,  for 
simplicity.  As  in  scalar  convection,  we  add  an  additional  term  to  control  the  shape  of  a 
short  pulse: 

d*0  =  c2d2J  +  d2xip 

or,  using  a  simple  time  discretization, 

-2f  +  fA  =  c2Atd2J  +  A td\yr  (3.3) 

It  was  seen  in  [SDHXLF]  that  the  addition  of  “Lattice  Confinement”  terms  in  the  form  of 
second  derivatives  of  a  fimction  that  decays  sufficiently  rapidly  away  from  the  centroid 
do  not  change  the  propagating  speed  (nor  the  total  amplitude)  of  a  propagating,  confined 
pulse.  The  same  is  true  here.  This  means,  of  course,  that  they  do  not  change  the  motion  of 
the  centroid  of  an  isolated  pulse. 

The  main  constraint  on  the  Confinement  term^/ ,  is  that  it  force  an  initial  isolated, 
propagating  compact  pulse  with  a  single  maximum  to  remain  compact  and  not  develop 
any  additional  maxima.  We  use: 

V"=^s;r-es;9’(r)  (3.4) 

In  this  term  O  has  the  form  given  by  Eqn.  (3.2)  in  terms  of  its  argument.  We  have 
defined 


Results  will  be  given  in  Sec.  4  using  this  form. 

An  important  feature  of  the  method  is  that  the  wave  do  not  suffer  a  “phase  shift” 
when  they  pass  through  each  other.  This  is  obvious  for  the  equation  we  want  to  simulate 
-  the  linear  wave  equation.  However,  the  Confinement  term  is  nonlinear.  Such  a  phase 
shift  would  show  up  as  a  kink  in  two  waves  in  2  or  3  dimensions  that  are  passing  through 
each  other,  and  can  be  studied  in  detail  in  1-D.  It  turns  out  that  there  is  no  kink,  to 
plottable  accuracy,  as  can  be  seen  in  the  plotted  scattering  results  in  Sec.  4.  Results  for 
the  centroid  trajectories  for  2  pulse  passing  through  each  other  in  1-D  are  presented  in 
Fig.  8.  There,  the  computed  centroids  are  plotted  as  solid  lines  and  the  exact  as  dashed 
(the  periodic  boundary  conditions  can  be  seen  in  the  former).  It  can  be  seen  that  there  is 
no  phase  shift  to  plottable  accuracy.  This  lack  of  nonlinear  interaction  persists,  according 
to  our  study,  in  the  limit  of  small  time  step  (2  orders  of  magnitude  smaller  that  that  of 
Fig.  8),  even  though  0(1 02)  confinement  corrections  were  applied  as  the  pulses  were 
overlapping.  We  attribute  this  to  the  existence  of  another  conserved  variable,  similar  to 
total  energy.  This  computation  was  done  by  Nick  Lynn  of  University  of  Tennessee  Space 
Institute  (UTSI)  under  separate  funding. 

Another  important  study  involves  the  pulse  speed  in  non-uniform  slowness  fields. 
This  problem  was  suggested  by  Alan  Newell.  For  an  initially  circular  pulse  in  2-D 
propagating  through  regions  where  slowness  varied  by  a  factor  of  5  across  the  circle, 
numerical  errors  in  the  speed  were  insignificant  to  plottable  accuracy.  This  was  verified 
by  varying  the  grid  size  by  a  factor  of  40  in  each  direction.  (This  was  done  by  Meng  Fan 
of  UTSI  under  separate  funding.) 


One  other  important  use  of  Lattice  Confinement  for  the  wave  equation  involves 
cases  with  multiple  grids  with  grid  interfaces.  These  are  used  in  the  missile  reflection 
problem  described  in  Sec.  4.  If  only  the  discretized  wave  equation  is  used,  with  no 
Confinement,  reflections  result  from  small  numerical  errors  at  the  grid  interfaces,  unless 
special  care  is  taken.  Lattice  Confinement  completely  overcomes  this  problem  and  the 
(single)  pulse  propagates  across  the  interface  with  no  reflections. 


Some  of  the  advantages  of  using  direct  difference  equations  instead  of  pde’s  for  a 
related  problem  can  be  found  in  [CSE].  There,  a  very  efficient  formulation  is  derived  for 
the  Helmholtz  equation  by  minimizing  the  L2  norm  of  the  error  of  waves  propagating 


with  fixed 


3.4  Vector  Potential  Formulation 

In  a  recent  promising  development  in  our  wave  equation  research,  we  have 
implemented  a  vector  potential,  A ,  as  the  basic  variable,  instead  of  a  scalar  representing 
the  magnitude  of  the  B  field.  Even  though  there  is  no  direct  effect  of  this  change  in 
classical  electromagnetics,  there  is  a  major  advantage  in  the  computation  of  pulse 
solutions.  This  is  because,  for  a  thin  pulse,  E  and  B  only  have  a  small  region  of  support, 
but  A  extends  throughout  the  field.  The  E  and  B  representation  of  the  pulse  is  a  spread 
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delta  function  across  the  width,  but  the  A  representation  is  a  step  function.  As  a  result,  it 
is  much  easier  to  capture  the  pulse  by  operating  on  A ,  since  it  is  topologically  “trapped” 
by  this  field.  The  argument  is  exactly  the  same  as  in  Vorticity  Confinement,  where  thin 
vortical  regions  ( a> )  are  trapped  in  a  velocity  field,  q ,  which  extends  throughout  space. 
(This  is  also  related  to  the  trapping  of  defects  in  other  field  theories).  The  resulting 
confinement  terms  are  exactly  the  same  in  the  two  cases,  with  A  corresponding  to  q  and 
B  corresponding  to  0) . 


The  formulation  is 


B  =  Vx  A 


S?A  =  c2V2DA  +  juS ,V2dA  +  sS ,V d  xb 

where  b  is  a  local  harmonic  mean  of  B  at  each  grid  point: 


yiAi 
r  n 


(3.5) 


and  where  8t  and  VD  denote  discrete  operators  and  Eqn.  3.5  was  used  with  a  sum  over 
N  neighboring  points,  labeled  £ .  In  the  formulation  the  Coulomb  gauge  was  used  so  that 
a  scalar  potential  does  not  appear  in  the  equation  for  B .  Also,  V  •  A  =  0  is  then  enforced. 
This  is  also  directly  analogous  to  the  incompressibility  condition,  V  •  q  =  0 ,  which  is 
continuing  on  this  formulation  for  the  surface  reflection  conditions.  We  are  currently 
investigating  the  properties  of  this  formulation  under  reflection. 


A  simpler  scalar  formulation  where  the  pulse  is  concentrated  near  a  discontinuity 
of  the  scalar,  as  in  the  acoustic  equation,  is  also  described  in  Sec.  4. 


4.  NUMERICAL  RESULTS 
4.1  Level  Set  Formulation 

We  have  applied  our  method  to  the  cases  of  2  and  3  dimensions  with  reflection 
and  refraction,  also  including  curved  boundaries.  Figure  1  shows  an  expanding  circle  in  a 
medium  of  index  of  refraction  1  in  such  a  setting  using  our  algorithm.  As  time  increases, 
more  and  more  reflections  take  place,  leading  to  wavefronts  that  almost  fill  up  spatial 
space.  Note  our  approach  not  only  handles  the  multiple  reflections,  which  lead  to 
complicated  multivalued  solutions,  but  resolves  well  the  wavefronts  which  have  grown  in 
length.  Figure  2  shows  the  cases  of  reflecting  to  the  curved  boundaries.  For  the  cases  with 
both  reflection  and  refraction,  in  figure  3  and  4,  we  set  the  two  materials  with  different 
index  of  reflection.  As  time  goes  on,  as  soon  as  our  expanding  circle  hits  the  other 
material,  both  reflection  and  refraction  take  place  simultaneously.  Finally,  Figure  5,  6  and 
7  show  the  three  dimensional  cases  with  different  initial  settings  (sphere,  ellipsoid,  torus). 

4.2.  Lattice  Confinement  Formulation 
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In  all  of  the  cases,  Lattice  Confinement  was  added  to  the  linear  wave  equation 
with  simple,  second  order  central  discretization. 

4.1  Reflection  from  Planar  Surface 

In  the  first  demonstration,  a  uniform  Cartesian  grid  was  used  with  reflection  on 
the  walls  of  a  square  region,  in  2-D  as  well  as  3-D.  In  general,  it  is  well  known  that  a 
“tail”  will  develop  behind  a  2-D  wave  front,  while  the  front  remains  sharp  (if  it  is  initially 
sharp).  In  the  2-D  computation  with  Confinement  we  can  keep  a  sharp  pulse  and  suppress 
the  tail.  This  tail  is  smooth  and  could,  if  desired,  be  computed  using  standard  CFD  on  the 
grid.  The  goal,  however,  is  to  show  that  a  pulse  can  propagate  over  long  distances  with 
Lattice  Confinement.  The  same  long  distance  propagation  was  observed  in  3-D  where 
there  were  genuine  pulse  solutions. 

In  Fig.9  2-D  results  are  presented  for  a  128x128  cell  grid.  It  can  be  seen  that  there 
is  no  perceptible  diffusion  of  the  pulse,  even  after  many  reflections.  These  results  were 
also  presented  in  [FWS]. 

4.2  Scattering  from  Missile 

This  computation  involved  multiple  body-fitted  grids  used.  The  first  missile 
computation  involved  a  planar  scalar  wave  impinging  on  the  nose  of  a  missile.  The  plane 
of  the  wave  is  parallel  to  the  longitudinal  axis  of  the  missile.  The  magnitude  of  the  scalar 
is  presented  in  planes  parallel  (i.e.,  the  symmetry  plane)  and  perpendicular  to  the 
longitudinal  axis  of  the  missile,  and  on  the  surface  of  the  missile  body.  Figures  10  and  1 1 
depict  the  intensity  of  the  planar  wave  before,  during,  and  after  reflection  from  the 
missile  surface  in  the  nose  area  from  the  side  and  from  the  front,  respectively. 

The  second  computation  involved  the  aft-end  of  the  missile.  Figure  12  depicts  the 
same  sequence,  but  in  the  fin  area,  as  seen  from  the  rear.  The  third  computation  involved 
the  entire  missile.  Figure  13  depicts  a  pulse  reflecting  from  the  nose  area  at 45°,  in  the 
symmetry  plane,  as  seen  from  the  side. 

Figure  14  depicts  part  of  the  overlapping  computational  grid.  The  confinement 
method  has  been  generalized  to  curvilinear  grids  for  this  case.  Besides  keeping  the  pulse 
sharp,  Confinement  also  prevents  reflections  from  the  grid  interface  regions  and  thus 
represents  an  important  capability. 

All  of  the  missile  results  involve  a  pulse  confined  to  about  3  grid  cells  except  in 
the  fine  grid  region  very  near  the  surface.  There,  the  method  reverts  to  standard  CFD. 
Also,  it  can  be  seen  in  Fig.  12  that,  unlike  in  geometrical  optics,  there  is  a  diffracted 
component.  This  will  most  likely  represent  a  diffracting  pulse  corresponding  to  the  width 
of  the  computational  one.  We  should  be  able  to  make  corrections  to  the  diffracted 
intensity  to  simulate  pulses  of  other  widths.  This  would  result  in  a  very  general  method 
able  to  treat  diffraction  to  first  order  (in  this  short  pulse  limit).  We  are  currently  studying 
this  possibility. 
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4.3  Scattering  from  Aircraft 

Next,  results  of  confinement  are  presented  for  scattering  from  an  aircraft  shape. 
Here,  a  uniform  Cartesian  grid  was  used  with  a  treatment  of  the  surface  boundary 
conditions  involving  interpolation.  A  level  set  description  of  the  body  was  used,  which 
was  early  derived  from  a  surface  definition  “STL”  file.  This  is  shown  in  Fig.  15.  With 
these  boundary  conditions,  it  can  be  seen  that  scattering  from  very  complex  shapes  can 
easily  be  computed  since  body  fitted  grids  are  not  required.  In  Fig.  16,  contours  of  the 
scalar  magnitude  representing  the  electromagnetic  field  in  a  cross  plane  near  the  back  of 
the  wing  are  shown  for  three  different  times. 

The  next  results  involve  the  new  vector  potential  computation.  There,  a  single 
contour  can  be  used  to  represent  each  part  of  the  wave  since  the  potential  extends 
throughout  the  field  and  has  a  definite  range  of  values  within  the  pulse.  Preliminary 
results  are  depicted  in  Fig.  17  for  two  times.  Further  work  is  progressing  on  the  new 
technique. 

The  last  results  depict  a  scalar  version  of  the  vector  potential  method  showing 
scattering  in  the  cross  plane  depicted  in  Fig.  15.  A  contour  representing  the  centroid  of 
the  pulse  is  shown  in  Fig.  18,  as  a  dashed  line.  Other  contours  showing  the  range  of 
values  of  the  scalar  are  also  shown. 
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